#!/usr/bin/env python
# -*- coding: UTF-8 -*-
'''
@Project ：V2 
@File    ：T2.py
@IDE     ：PyCharm 
@Author  ：郭星
@Date    ：2025/10/23 16:04 
'''
# -*- coding: utf-8 -*-
"""
gurobi_gas_opt.py
---------------------------------
用 Gurobi 构建“标气生产车间（清洗→充装→摇匀→分析）”的 MILP 模型，
用于在小规模算例上求真最优（或强下界），以验证启发式/元启发式算法（如改进禁忌搜索）的有效性。

⚠️ 说明：
1) 本模型支持：
   - 单/多台充装设备（配气机）。
   - 单台清洗设备、单台摇匀设备（如需多台，可按充装设备的建模方式扩展）。
   - 分析阶段“分组并测 + 同步起止 + 组时长 = 组内各组分分析时长的最大值”。
   - 分析设备能力矩阵 can[d][c]（例如：设备1能AB、设备2能BC、设备3能CD）。
2) 目标：最小化最大完工时间（makespan）。
3) 为了使模型规模可控，建议小规模验证：瓶数 ≤ 6、组分 ≤ 5、组上限 ≤ 3、设备数 ≤ 5。

依赖：gurobipy
命令行：python gurobi_gas_opt.py
"""

from typing import List, Dict, Tuple
import itertools
import math

try:
    import gurobipy as gp
    from gurobipy import GRB
except Exception as e:
    raise RuntimeError("请先安装 Gurobi 并配置许可（pip install gurobipy）。当前错误：%s" % str(e))


# ============ 辅助：构造一个可跑的小型示例（你也可以替换为自己的真实数据） ============
def build_toy_instance():
    """
    返回一个字典，包含模型需要的全部数据。
    设定：
    - 2 瓶产品：瓶1需要 ABC，瓶2需要 BCD；底气 N2 作为最后一个组分。
    - 充装设备 1 台（你可改为 >1）。
    - 分析设备 3 台：D1 能 AB，D2 能 BC，D3 能 CD。
    - 清洗、充装、摇匀、分析的时间单位统一为“分钟”。
    """
    instance = {}

    # 瓶集合与组分集合（最后一个是底气 N2）
    B = 2                    # 瓶数
    components = ["A", "B", "C", "D", "N2"]
    C = len(components)
    bottom = C - 1           # 约定：最后一个组分是底气

    # 每瓶需要的有效组分（不含底气），可用于校验或限制（本模型不强制使用该列表，按时间>0视为需要）
    need_mask = [[1,1,1,0,1],   # 瓶1：A,B,C + N2（最后一个1表示需要底气）
                 [0,1,1,1,1]]   # 瓶2：  B,C,D + N2

    # 清洗 / 充装 / 摇匀 / 分析时间
    clean_time = [5, 5]                              # 每瓶清洗时长
    shake_time = [3, 3]                              # 每瓶摇匀时长
    # 充装时间：矩阵 [B][C]
    # 若某瓶某组分不需要，时间也可以填较小值（但依然会被建模）。
    fill_time = [
        [8, 6, 7, 0, 4],  # 瓶1：A,B,C,D,N2（D=0 意味着不需要该组分，仍会顺序约束避免N2在前）
        [0, 7, 6, 8, 4],  # 瓶2：A,B,C,D,N2
    ]

    # 分析设备能力与时间：
    # 3 台分析设备：
    #  - D1 能 A,B
    #  - D2 能 B,C
    #  - D3 能 C,D
    D = 3
    can = [
        [1, 1, 0, 0, 0],  # D1: A,B
        [0, 1, 1, 0, 0],  # D2: B,C
        [0, 0, 1, 1, 0],  # D3: C,D
    ]
    # 对应的分析时间（若 can[d][c]==0 则该时间无效）
    anal_time = [
        [3, 2, 0, 0, 0],  # D1：A=3, B=2
        [0, 2, 1, 0, 0],  # D2：B=2, C=1
        [0, 0, 2, 3, 0],  # D3：C=2, D=3
    ]

    # 设备数量：
    F = 1  # 充装设备数量（可改 1→2+）
    # 清洗/摇匀：这里用单机，若你要多机，可仿照充装设备进行“指派+两两不重叠”建模。
    CL = 1  # 清洗设备台数
    SK = 1  # 摇匀设备台数（目前按单机建模，如需多机自行扩展）

    # 分析分组上限（每瓶最多分几组做分析，同时一瓶同一时刻只能接一套分析管路）
    Gmax = 3

    instance.update(dict(
        B=B, C=C, D=D, F=F, CL=CL, SK=SK, Gmax=Gmax,
        components=components, bottom=bottom,
        clean_time=clean_time, fill_time=fill_time, shake_time=shake_time,
        can=can, anal_time=anal_time, need_mask=need_mask
    ))
    return instance


# ============ 主建模函数 ============
def build_and_solve(inst: Dict, time_limit: int = 60, mip_gap: float = 1e-4):
    """
    根据给定实例构建 Gurobi 模型并求解。
    :param inst: 由 build_toy_instance() 或你自己的数据结构返回
    :param time_limit: Gurobi 求解时间上限（秒）
    :param mip_gap: MIPGap 收敛阈值
    :return: (model, solution_dict)
    """
    B = inst['B']
    C = inst['C']
    Dv = inst['D']
    F = inst['F']
    CL = inst['CL']
    SK = inst['SK']
    Gmax = inst['Gmax']
    components = inst['components']
    bottom = inst['bottom']

    clean_time: List[float] = inst['clean_time']
    fill_time: List[List[float]] = inst['fill_time']
    shake_time: List[float] = inst['shake_time']
    can: List[List[int]] = inst['can']
    anal_time: List[List[float]] = inst['anal_time']

    # ========== 计算一个足够大的 M（上界） ==========
    T_clean = sum(clean_time)
    T_fill = sum(sum(max(ft, 0.0) for ft in row) for row in fill_time)
    T_shake = sum(shake_time)
    # 分析时间的保守上界：对每瓶，最多把所有有效组分串行做一次“最大时长”，乘设备数上界不必要，这里简化。
    T_anal = 0.0
    for b in range(B):
        # 取各组分的 max_d anal_time[d][c]
        for c in range(C-1):  # 不含底气
            best = max(anal_time[d][c] if can[d][c] else 0.0 for d in range(Dv))
            T_anal += best
    M = 10.0 * (T_clean + T_fill + T_shake + T_anal + 1.0)  # 乘个系数，防止边界

    # ========== 模型 ==========
    m = gp.Model("GasMixShop_MILP")

    # ---------- 变量：时间 ----------
    S_clean = m.addVars(B, vtype=GRB.CONTINUOUS, lb=0.0, name="S_clean")
    E_clean = m.addVars(B, vtype=GRB.CONTINUOUS, lb=0.0, name="E_clean")

    S_fill = m.addVars(B, C, vtype=GRB.CONTINUOUS, lb=0.0, name="S_fill")
    E_fill = m.addVars(B, C, vtype=GRB.CONTINUOUS, lb=0.0, name="E_fill")

    S_shake = m.addVars(B, vtype=GRB.CONTINUOUS, lb=0.0, name="S_shake")
    E_shake = m.addVars(B, vtype=GRB.CONTINUOUS, lb=0.0, name="E_shake")

    # 每瓶最后一个充装完成的时刻（用于约束摇匀开始）
    E_fill_last = m.addVars(B, vtype=GRB.CONTINUOUS, lb=0.0, name="E_fill_last")

    # 分析分组：每瓶至多 Gmax 组
    S_grp = m.addVars(B, Gmax, vtype=GRB.CONTINUOUS, lb=0.0, name="S_grp")
    E_grp = m.addVars(B, Gmax, vtype=GRB.CONTINUOUS, lb=0.0, name="E_grp")
    T_grp = m.addVars(B, Gmax, vtype=GRB.CONTINUOUS, lb=0.0, name="T_grp")

    # 完工时间 & 目标 makespan
    C_b = m.addVars(B, vtype=GRB.CONTINUOUS, lb=0.0, name="C_b")
    MKS = m.addVar(vtype=GRB.CONTINUOUS, lb=0.0, name="Makespan")

    # ---------- 变量：设备指派与不重叠 ----------
    # 清洗/摇匀：这里按单台设备处理 → 用两两不重叠的 0/1 排序变量
    o_clean = m.addVars(B, B, vtype=GRB.BINARY, name="ord_clean")   # o_clean[b,b']=1 表示 b 在 b' 之前清洗
    o_shake = m.addVars(B, B, vtype=GRB.BINARY, name="ord_shake")

    # 充装设备指派
    a_fill = m.addVars(B, C, F, vtype=GRB.BINARY, name="a_fill")    # a_fill[b,c,f]=1 表示(b,c)在设备f上充装
    # 同一设备上的不重叠
    u_fill = m.addVars(B, C, B, C, F, vtype=GRB.BINARY, name="u_fill")  # u_fill[i,j,f]=1 表示 i 在 j 之前（在同一台 f 上）

    # 同瓶充装“串行”约束（即使多机也不允许一瓶多组分并行）
    k_serial = m.addVars(B, C, C, vtype=GRB.BINARY, name="serial_fill")  # k_serial[b,c,c']=1 表示 c 在 c' 之前

    # 分析阶段：分组指示 / 组分分配 / 设备使用
    v_grp = m.addVars(B, Gmax, vtype=GRB.BINARY, name="v_grp")           # 该组是否存在（含至少1个组分）
    z = m.addVars(B, C-1, Gmax, vtype=GRB.BINARY, name="z")              # 组分 c(不含N2) 分配到组 g
    x = m.addVars(B, C-1, Gmax, Dv, vtype=GRB.BINARY, name="x")          # 若 z=1，则必须选择且仅选1台可测设备
    y = m.addVars(B, Gmax, Dv, vtype=GRB.BINARY, name="y")               # 组 g 是否占用设备 d（存在某个 c 用了该 d）

    # 设备 d 在不同组间的不重叠（只在 y=1 的组对之间生效）
    o_dev = m.addVars(B, Gmax, B, Gmax, Dv, vtype=GRB.BINARY, name="ord_dev")

    # 同瓶不同组不重叠（同一瓶同一时刻只能接一套管路）
    o_grp = m.addVars(B, Gmax, Gmax, vtype=GRB.BINARY, name="ord_grp")

    # ========== 约束：清洗阶段 ==========
    for b in range(B):
        m.addConstr(E_clean[b] == S_clean[b] + clean_time[b], name=f"C_clean_dur[{b}]")

    # 单机不重叠（清洗）
    for b in range(B):
        for bp in range(B):
            if b == bp:
                continue
            m.addConstr(S_clean[bp] >= E_clean[b] - M * (1 - o_clean[b, bp]), name=f"C_clean_seq1[{b},{bp}]")
            m.addConstr(S_clean[b]  >= E_clean[bp] - M * (o_clean[b, bp]),  name=f"C_clean_seq2[{b},{bp}]")

    # ========== 约束：充装阶段 ==========
    # 时长
    for b in range(B):
        for c in range(C):
            m.addConstr(E_fill[b, c] == S_fill[b, c] + fill_time[b][c], name=f"C_fill_dur[{b},{c}]")
            # 充装不得早于清洗完成
            m.addConstr(S_fill[b, c] >= E_clean[b], name=f"C_fill_after_clean[{b},{c}]")

    # 每个充装工序必须指派到且仅指派到一台充装设备
    for b in range(B):
        for c in range(C):
            m.addConstr(gp.quicksum(a_fill[b, c, f] for f in range(F)) == 1, name=f"C_fill_assign[{b},{c}]")

    # 同一设备上的不重叠（条件：两工序都在该设备）
    for f in range(F):
        for b, c, bp, cp in itertools.product(range(B), range(C), range(B), range(C)):
            if (b, c) >= (bp, cp):
                continue
            # (b,c) 在 (bp,cp) 之前
            m.addConstr(S_fill[bp, cp] >= E_fill[b, c] - M * (1 - u_fill[b, c, bp, cp, f])
                        - M * (1 - a_fill[b, c, f]) - M * (1 - a_fill[b, cp, f]),
                        name=f"C_fill_noovl1[{b},{c},{bp},{cp},{f}]")
            # (bp,cp) 在 (b,c) 之前
            m.addConstr(S_fill[b, c] >= E_fill[bp, cp] - M * (u_fill[b, c, bp, cp, f])
                        - M * (1 - a_fill[b, c, f]) - M * (1 - a_fill[bp, cp, f]),
                        name=f"C_fill_noovl2[{b},{c},{bp},{cp},{f}]")

    # 同瓶充装串行（不允许并行）+ 底气必须在该瓶最后
    for b in range(B):
        for c in range(C):
            for cp in range(C):
                if c == cp:
                    continue
                # 同瓶串行：二选一
                m.addConstr(S_fill[b, cp] >= E_fill[b, c] - M * (1 - k_serial[b, c, cp]),
                            name=f"C_fill_serial1[{b},{c},{cp}]")
                m.addConstr(S_fill[b, c]  >= E_fill[b, cp] - M * (k_serial[b, c, cp]),
                            name=f"C_fill_serial2[{b},{c},{cp}]")

        # 底气最后：对每个非底气组分 c，强制 E_fill[b,c] <= S_fill[b,bottom]
        for c in range(C-1):
            m.addConstr(E_fill[b, c] <= S_fill[b, bottom], name=f"C_fill_bottom_last[{b},{c}]")

    # 计算每瓶最后一个充装完工时刻，用于约束摇匀开始
    for b in range(B):
        for c in range(C):
            m.addConstr(E_fill_last[b] >= E_fill[b, c], name=f"C_fill_lastmax[{b},{c}]")

    # ========== 约束：摇匀阶段（单机） ==========
    for b in range(B):
        m.addConstr(E_shake[b] == S_shake[b] + shake_time[b], name=f"C_shake_dur[{b}]")
        m.addConstr(S_shake[b] >= E_fill_last[b], name=f"C_shake_after_fill[{b}]")

    for b in range(B):
        for bp in range(B):
            if b == bp:
                continue
            m.addConstr(S_shake[bp] >= E_shake[b] - M * (1 - o_shake[b, bp]), name=f"C_shake_noovl1[{b},{bp}]")
            m.addConstr(S_shake[b]  >= E_shake[bp] - M * (o_shake[b, bp]),  name=f"C_shake_noovl2[{b},{bp}]")

    # ========== 约束：分析阶段（分组并测 + 同步起止 + 设备能力） ==========
    # 组存在与组分覆盖
    for b in range(B):
        for c in range(C-1):  # 不含底气
            m.addConstr(gp.quicksum(z[b, c, g] for g in range(Gmax)) == 1, name=f"C_group_cover[{b},{c}]")
        # 组存在指示：若某组有任何组分，则 v_grp=1；否则 v_grp=0
        for g in range(Gmax):
            for c in range(C-1):
                m.addConstr(v_grp[b, g] >= z[b, c, g], name=f"C_group_exist_lb[{b},{g},{c}]")
            m.addConstr(v_grp[b, g] <= gp.quicksum(z[b, c, g] for c in range(C-1)), name=f"C_group_exist_ub[{b},{g}]")

    # 每个(瓶,组,组分)必须选择**且仅选择**一台能测该组分的设备
    for b in range(B):
        for c in range(C-1):
            for g in range(Gmax):
                # 只能从 can[d][c]==1 的设备中选
                m.addConstr(gp.quicksum(x[b, c, g, d] for d in range(Dv) if can[d][c] == 1) == z[b, c, g],
                            name=f"C_group_device_choose[{b},{c},{g}]")
                # 对 can==0 的设备强制 x=0（可省略，因为上式已经只在 can==1 的 d 上求和，但写清晰）
                for d in range(Dv):
                    if can[d][c] == 0:
                        m.addConstr(x[b, c, g, d] == 0, name=f"C_forbid_incap[{b},{c},{g},{d}]")

    # 每组同一台设备同时最多测 1 个组分（避免“一机多活”）
    for b in range(B):
        for g in range(Gmax):
            for d in range(Dv):
                m.addConstr(gp.quicksum(x[b, c, g, d] for c in range(C-1)) <= 1, name=f"C_one_comp_per_dev[{b},{g},{d}]")
                # 定义 y = 是否使用设备 d（恰为上述求和）
                m.addConstr(y[b, g, d] == gp.quicksum(x[b, c, g, d] for c in range(C-1)), name=f"C_y_def[{b},{g},{d}]")

    # 组持续时间 = 组内被选设备的分析时长的“最大值” → 线性化为一组下界
    for b in range(B):
        for g in range(Gmax):
            # 若组不存在，则 T_grp=0（也可仅上界限制，目标会自动压到0；这里加上更严格些）
            m.addConstr(T_grp[b, g] <= M * v_grp[b, g], name=f"C_Tgrp_zero_if_empty[{b},{g}]")
            for c in range(C-1):
                for d in range(Dv):
                    if can[d][c] == 1:
                        # T_grp >= anal_time[d][c] * x[b,c,g,d]
                        m.addConstr(T_grp[b, g] >= anal_time[d][c] * x[b, c, g, d],
                                    name=f"C_Tgrp_maxlb[{b},{g},{c},{d}]")
            # 起止
            m.addConstr(E_grp[b, g] == S_grp[b, g] + T_grp[b, g], name=f"C_grp_dur[{b},{g}]")
            # 不能早于摇匀完成（若该组存在）
            m.addConstr(S_grp[b, g] >= E_shake[b] - M * (1 - v_grp[b, g]), name=f"C_grp_after_shake[{b},{g}]")

    # 同瓶不同组不重叠（同一瓶同一时刻只能接一套分析管路）
    for b in range(B):
        for g in range(Gmax):
            for gp2 in range(Gmax):
                if g == gp2:
                    continue
                m.addConstr(S_grp[b, gp2] >= E_grp[b, g] - M * (1 - o_grp[b, g, gp2])
                            - M * (1 - v_grp[b, g]) - M * (1 - v_grp[b, gp2]),
                            name=f"C_grp_noovl1[{b},{g},{gp2}]")
                m.addConstr(S_grp[b, g] >= E_grp[b, gp2] - M * (o_grp[b, g, gp2])
                            - M * (1 - v_grp[b, g]) - M * (1 - v_grp[b, gp2]),
                            name=f"C_grp_noovl2[{b},{g},{gp2}]")

    # 同一设备在不同(瓶,组)间不重叠（条件：两个组都使用该设备）
    for d in range(Dv):
        for b, g, bp, gp2 in itertools.product(range(B), range(Gmax), range(B), range(Gmax)):
            if (b, g) >= (bp, gp2):
                continue
            m.addConstr(S_grp[bp, gp2] >= E_grp[b, g] - M * (1 - o_dev[b, g, bp, gp2, d])
                        - M * (1 - y[b, g, d]) - M * (1 - y[bp, gp2, d]),
                        name=f"C_dev_noovl1[{b},{g},{bp},{gp2},{d}]")
            m.addConstr(S_grp[b, g] >= E_grp[bp, gp2] - M * (o_dev[b, g, bp, gp2, d])
                        - M * (1 - y[b, g, d]) - M * (1 - y[bp, gp2, d]),
                        name=f"C_dev_noovl2[{b},{g},{bp},{gp2},{d}]")

    # ========== 完工与目标 ==========
    for b in range(B):
        # 每瓶的完工时间 C_b ≥ 所有组的结束时间；也 ≥ 摇匀结束（以防某瓶异常没有组）
        m.addConstr(C_b[b] >= E_shake[b], name=f"C_cb_after_shake[{b}]")
        for g in range(Gmax):
            m.addConstr(C_b[b] >= E_grp[b, g], name=f"C_cb_after_grp[{b},{g}]")
        m.addConstr(MKS >= C_b[b], name=f"C_makespan_def[{b}]")

    m.setObjective(MKS, GRB.MINIMIZE)

    # ========== 参数 & 求解 ==========
    m.Params.TimeLimit = time_limit
    m.Params.MIPGap = mip_gap
    # m.Params.OutputFlag = 0  # 如需静默求解可打开

    m.optimize()

    # ========== 提取解 ==========
    sol = {}
    if m.status in [GRB.OPTIMAL, GRB.TIME_LIMIT, GRB.SUBOPTIMAL]:
        sol['obj'] = m.objVal if m.SolCount > 0 else None
        sol['status'] = m.status
        sol['S_clean'] = {b: S_clean[b].X for b in range(B)}
        sol['E_clean'] = {b: E_clean[b].X for b in range(B)}
        sol['S_fill'] = {(b, c): S_fill[b, c].X for b in range(B) for c in range(C)}
        sol['E_fill'] = {(b, c): E_fill[b, c].X for b in range(B) for c in range(C)}
        sol['S_shake'] = {b: S_shake[b].X for b in range(B)}
        sol['E_shake'] = {b: E_shake[b].X for b in range(B)}
        sol['S_grp'] = {(b, g): S_grp[b, g].X for b in range(B) for g in range(Gmax)}
        sol['E_grp'] = {(b, g): E_grp[b, g].X for b in range(B) for g in range(Gmax)}
        sol['T_grp'] = {(b, g): T_grp[b, g].X for b in range(B) for g in range(Gmax)}
        sol['C_b'] = {b: C_b[b].X for b in range(B)}
        sol['MKS'] = MKS.X

        # 分析组的结构
        sol['z'] = {(b, c, g): int(round(z[b, c, g].X)) for b in range(B) for c in range(C-1) for g in range(Gmax)}
        sol['x'] = {(b, c, g, d): int(round(x[b, c, g, d].X)) for b in range(B) for c in range(C-1) for g in range(Gmax) for d in range(Dv)}
        sol['y'] = {(b, g, d): int(round(y[b, g, d].X)) for b in range(B) for g in range(Gmax) for d in range(Dv)}

    else:
        sol['status'] = m.status
        sol['obj'] = None

    return m, sol


# ============ 打印可读日程（便于对比启发式结果） ============
def pretty_print_solution(inst: Dict, sol: Dict):
    B, C, Dv, Gmax = inst['B'], inst['C'], inst['D'], inst['Gmax']
    components, bottom = inst['components'], inst['bottom']
    anal_time, can = inst['anal_time'], inst['can']

    print("\n================= 求解结果 =================")
    print(f"状态: {sol['status']}  目标(Makespan): {sol.get('MKS')}")
    print("--------------------------------------------")
    for b in range(B):
        print(f"[瓶 {b}] 清洗: {sol['S_clean'][b]:.1f} → {sol['E_clean'][b]:.1f}")
        # 充装（按开始时间排序输出）
        pairs = [((b, c), sol['S_fill'][b, c], sol['E_fill'][b, c]) for c in range(C)]
        pairs.sort(key=lambda t: t[1])
        for (bb, c), s, e in pairs:
            tag = " (底气)" if c == bottom else ""
            print(f"    充装 组分 {components[c]}{tag}: {s:.1f} → {e:.1f}")
        print(f"    摇匀: {sol['S_shake'][b]:.1f} → {sol['E_shake'][b]:.1f}")
        # 分析组
        grp_list = [((b, g), sol['S_grp'][b, g], sol['E_grp'][b, g], sol['T_grp'][b, g]) for g in range(Gmax)]
        grp_list = [t for t in grp_list if t[1] + t[2] + t[3] > 1e-6]  # 过滤空组
        grp_list.sort(key=lambda t: t[1])
        for (bb, g), s, e, T in grp_list:
            print(f"    分析 组 g={g}: {s:.1f} → {e:.1f}  (T={T:.1f})")
            # 列出组内各组分及所用设备
            for c in range(C-1):
                for d in range(Dv):
                    if sol['x'][b, c, g, d] == 1:
                        print(f"        组分 {components[c]} → 设备 d={d}  (t={anal_time[d][c]}; 可用={can[d][c]})")
        print(f"    完工 C_b={sol['C_b'][b]:.1f}")
        print("--------------------------------------------")
    print(f"总完工时间 Makespan = {sol.get('MKS')}")
    print("============================================\n")


# ============ 主程序入口 ============
def main():
    # 你也可以把 build_toy_instance() 换成读取 JSON/CSV 的自定义函数
    inst = build_toy_instance()
    model, sol = build_and_solve(inst, time_limit=60, mip_gap=1e-4)

    if sol.get('obj') is not None:
        pretty_print_solution(inst, sol)
    else:
        print("未获得可用解，状态码：", sol['status'])


if __name__ == "__main__":
    main()
